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We study the stability of a Bose-Einstein condensate of harmonically 
trapped atoms with negative scattering length, specifically ''Li. Our method 
is to solve the time-dependent nonlinear Schrodinger equation numerically. 



For an isolated condensate, with no gain or loss, we find that the system is 
stable (apart from quantum tunneling) if the particle number N is less than a 
critical number Nc- For N > Nc, the system collapses to high-density clumps 
in a region near the center of the trap. The time for the onset of collapse is 



< 
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■ on the order of 1 trap period. Within numerical uncertainty, the results are 

■ consistent with the formation of a "black hole" of infinite density fluctuations, 
00 ' s-s predicted by Ueda and Huang We obtain numerically Nc ~ 1251. We 

I then include gain-loss mechanisms, i.e., the gain of atoms from a surrounding 

I "thermal cloud", and the loss due to two- and three-body collisions. The 

' number N now oscillates in a steady state, with a period of about 145 trap 

g ■ periods. We obtain Nc ~ 1260 as the maximum value in the oscillations. 

C ; I. INTRODUCTION AND SUMMARY 

o 
o 

Bose-Einstein condensation has been observed in magnetically trapped dilute vapors of 
the alkah elements ^'''Rb [|l|, ^^Na 0, ''Li and 0. At the nanodegree temperatures 
of these experiments, the systems would have frozen solid long ago were they in free space. 
5^ I In the confining trap, however, zero-point motion keeps the atoms apart, and the systems 
remain gaseous. The case of ^Li is special, however, in that the interatomic interaction is 
predominantly attractive, as indicated by a negative scattering length. Thus, the condensate 
in ^Li should be less stable than the other cases. The purpose of this paper is to study the 
nature of the instability, its onset, and manifestations. 

The object of study is the condensate wave function ilj{r,t), which gives the proba- 
bility amplitude for annihilating one particle in the condensate at r at time t. We use the 
time-dependent Gross-Pitaevskii (GP) equation, or nonlinear Schrodinger equation (NLSE), 
which corresponds to a mean-field approximation: 
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where a is a negative scattering length, and the external potential is taken to be harmonic: 



V{r) = ^muj^r^ (2) 

This defines a characteristic length do, the width of the unperturbed ground-state wave 
function: 



do = ^— (3) 

V muj 

The number of condensate particles enters through the normalization 

N = [ rfVl^p (4) 



which is a constant of the motion. Another constant of the motion is the Hamiltonian 



H = d^r 



fj2 TJ 



(5) 



The parameters used in ^Li experiments correspond to [Q 



a = —1.45 nm 
do ^ 3.16 fim 

CO fti 908 (6) 

where u is taken to be approximately equal to the geometric mean of the three circular 
frequencies of the experimental trap. Equation (|l|) does not take into account the coupling 
between the condensate and the "thermal cloud" of uncondensed atoms, nor does it describe 
the loss of atoms from the trap due to collisions. These effects will be considered later. 

In free space, the NLSE with attractive interactions has an instability known in plasma 
physics as "self-focusing" |^ , whereby the initial wave function develops a singularity in finite 
time, corresponding to a local collapse to a state of infinite density. In an external trap, 
however, the ground state is apparently stable, as long as is not too large, as indicated by 
variational calculations using a Gaussian trial wave function 0. The width of the Gaussian 
narrows with increasing N, and collapses to zero at some critical value A''^. This conclusion 
is borne out by other studies If^^TTI . In particular, Kim and Zubarev obtain an exact upper 
bound for N^, which for a Gaussian wave function gives 

N, ^ 0.671 ^ (7) 

\a\ 

Ueda and Leggett obtain the upper bound in a variational calculation. For the ^Li pa- 
rameters (6), this formula gives Nc ^ 1463. In the numerical calculations described later, 
we obtain 

N, = 0.574-^ (8) 
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or Nc = 1251 for the ^Li parameters. This number becomes shghtly larger when gain and 
loss effects are taken into account. 

Actually, even ignoring collisional loss, the condensate is only metastable for N < N^, 
for it can decay via quantum tunneling. This effect is not described by the NLSE, al- 
though an approximate decay amplitude can be obtained by continuation of the NLSE to 
imaginary time. Kagan et al. estimated the amplitude by calculating an overlap inte- 
gral between initial and final states represented by Gaussian wave functions, and obtained 
(df/cii)^^/^, where di and df are respectively the widths of the initial and final wave functions. 
Through numerical calculations, Shuryak found that the decay rate is proportional to 
exp[— const. (Ai'c — A^)]. Stoof [|T^ wrote down a WKB formula for the decay rate, but did 
not explicitly evaluate it. Using variational wave functions, Ueda and Leggett obtained 
a rate proportional to exp[— const. (Ai'c — Ny^^]. These estimates indicate that the tunneling 
probability is negligible unless N ^ N^. However, the WKB approximation breaks down 
in this neighborhood, and reliable calculations become difficult. In experimental terms, it 
is also difficult to observe tunneling in the present system, because the effect is masked 
by collisional loss, especially near N = Nc [|1^]. For these reasons, we shall not calculate 
the tunneling amplitude in this paper; but we will give a qualitative discussion of the 
phenomenon, in as far as it pertains to the instability of the system. 

Ueda and Huang formulated an approach to the stability problem, including tun- 
neling, in terms of the Feynman path integral for a transition amplitude. In a saddle-point 
approximation to the path integral, one obtains the NLSE as the saddle-point condition, 
whose continuation into imaginary time yields the tunneling amplitude. Starting with a 
Gaussian initial wave function, they assume that it remains Gaussian, in order to analyze 
the problem analytically. In this "Gaussian approximation", it was found that, as soon as 

> Ac, the system begins to collapse locally, in a region at the center of the trap. The 
size of this region grows as A^ increases. By putting the NLSE in hydrodynamic form, one 
can see that the pressure becomes negative inside this "black hole" , and, consequently, there 
would be infinite density fluctuations. 

In the present work, we remove the Gaussian approximation by performing numerical 
calculations. We verify that the picture presented by Ueda and Huang is correct. For A^ > 
A'c , a "black-hole" does appear at the center of the trap, within which the density has sharp 
local maxima, whose heights grow as the grid size for the numerical computation is decreased. 
This is consistent with the expectation that the density becomes divergent in the continuum 
limit. Finally, we take into account gain and loss mechanisms by adding phenomenological 
terms to the NLSE. We are able to study in some detail the growth-collapse cycles of the 
condensate, as anticipated in earlier work [|T^,|T^. 



II. TUNNELING AND COLLAPSE 

To give an overview of the mechanisms for instability, we think of the condensate wave 
function at each spatial point as a "coordinate," which moves in an effective potential, 
consisting of the last two terms in the Hamiltonian (|^), plus fluctuation corrections arising 
from the term tp*'V'^ip. The system is classically stable when this motion is bounded at all 
spatial points. We immediately see that the attractive interaction term —^{tp*tpy can lead 
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to instability, for it decreases without bound when the density increases. The question is 
whether there is an energy barrier guarding against the free fall. 

A precise definition for the effective potential can be modeled after that for a spatially 



uniform system in quantum field theory [|19|. For a qualitative description, however, it is 



much simpler to use the "Gaussian approximation" introduced by Ueda and Huang |]T6 
where we assume that the wave function is Gaussian, given by its initial form 



^o(r) =Coexp(-rVci2) 

where 

Co = n-'/^d-'/'VN 
We introduce a width parameter a by putting 

d = - 

'a 



(9) 



(10) 



(11) 



Thus, a = 1 corresponds to the ground-state eigenfunction of the harmonic oscillator. With 
this, we have 



^2 

2m 



where 



2mdl 



3a - I — 

do 



a 



(12) 



(13) 



and in this approximation the effective potential is 



nrW = [V{r) + x{r)]\ij\'- 



Uo 



(14) 



This is depicted on the right panel of Fig. |1| for different r. There is an energy barrier with 
height 



Wr = ^[V{r)+x{r)f 



(15) 



which is nonzero at r = 0, because x(0) 7^ 0. At large the effective potential tends to 
—00. In reality, of course, the NLSE ceases to be valid somewhere along the drop, for other 
physical effects, such as solidification, come in. 

A typical initial wave function ^/'o(^) is shown on the left panel of Fig. |l|. At a given r, 
we can measure the wave function on the left panel, and transfer it to the horizontal axis on 
the right panel. If it lands to the left of the barrier maximum for the particular r, then the 
system at that point is classically stable, but can decay via quantum tunneling as indicated 
by the classically forbidden path A —>■ B. If it lands to the right, the system at that r will 
rapidly collapse to a state of high density. 
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Since ^jjoii") has a maximum at r = with value proportional to \/N, the system may 
only decay via tunneling at any r, if it is classically stable at r = 0. Otherwise, the system 
in a region about r = will rapidly collapse. The condition for stability against collapse is 
therefore 

QoiMO)) < Wo (16) 

which corresponds to < Nc, with 

iVe = -,/^^ = ^^ (17) 

8 V « v« \0'\ 

The time for the onset of local collapse is expected to be greater than the time needed to 
establish a quasi- stationary initial wave function. Thus, a should correspond to the best 
Gaussian approximation to that wave function. This expectation is verified by numerical 
calculations described below, which also show that the Gaussian approximation is very good, 
up to the onset of collapse. Of course, a ceases to have meaning once the collapse begins, 
for > Nc- The numerical calculations for = 1248, which is just below N^., give 

a = 1.79 (18) 

As mentioned earlier, we will not calculate the tunneling amplitude in this paper, but 
concentrate on numerical solutions of the NLSE in real time. 



III. STATIONARY STATE 



We seek spherically symmetric solutions to the NLSE (Q), for the parameters (6) for the 
^Li experiments. For convenience we put 



,/ N IN u(r,t) 
with the boundary condition u{0) = 0, and the normalization condition 

1 r°° 

— / dr\u{r,t)\^ = 1 (20) 

From now on, we shall measure distance in units of do, and time in units of 2c(j^^. To do 
this without introducing new symbols for dimensionless distance and time, we formally put 
do = uj/2 = l. The NLSE (g) then takes the form 



where 



.du { 2 2G|mP . 



G^^^JL (22) 
do 2180 ^ ' 
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The initial condition with Gaussian form is 

u{r,0) = Morexp ^"2*^^^^ ^^^^ 

where 

uo = 27r-i/^«3/4 (24) 

We calculate u numerically as a complex function, using an implicit difference equation 
algorithm due to Goldberg et al. pO| . 



To find out how various initial states relax to a stationary state, we first solve the 
equation with an initial wave function uniform within a certain radius, and zero outside. 
For = 1142, which is somewhat below the critical number, the time development is shown 
in Fig. 0(a). We see that the modulus r|'i/'| appears to be Gaussian-like, except for small 
local fluctuations, after only about one-hundredth of an oscillator cycle. In Fig. ^(b), we 
start with a Gaussian with a = 1, with the same A^. This corresponds to the stationary 
wave function in the absence of interatomic interaction. In a few oscillator cycles, the width 
narrows to a stationary value. If we initially choose a near the stationary value, the location 
of the peak of rl-?/;! will slowly oscillate about some mean position. Fig. ^(c) shows the case 
when the ratio of the standard deviation to the mean of the peak location is minimized. 
This corresponds to a = 1.48, for this choice of N. The wave function in this case is sensibly 
stationary, and the Gaussian approximation is good. 

In Fig. ^, we plot the ratio of the standard deviation to the mean of the peak location for 
N = 1248 (just below critical) against a, and exhibit the minimum at a = 1.79, as quoted 
earlier in (p!8|). For a quantitative measure of how good the Gaussian approximation is, we 
calculate the absolute value of the overlap between u{r,t) and u{r,0) and find that it lies 
between 1 and 0.9795, with an average of 0.9906, for t between and 50 {t in units of 2uj~^). 



IV. LOCAL COLLAPSE 



To find Nc, we examine the time development for different A^, to look for the onset of local 
collapse. The signature of the collapse is the occurrence of a minimum in the wave function 
within a small distance from the origin as shown in Fig. | for a = 1.79 and = 1308. The 
initial wave function is stationary for a while, but begins to narrow after t ~ 0.2 (Fig. ^(a)), 
and a dip occurs at t ~ 0.7 (Fig. §(b)), at r ^ 2.3, which describes the implosion of a 
shell of particles at that radius. The first local minimum moves inwards with time and it 
eventually gets very close to r = 0, at t ~ 0.9, as shown in Fig. §(c). The density in the 
inner region can grow only by taking particles from the outside, since the number of particles 
is conserved. An examination of the phase of the wave function, whose gradient gives the 
superfiuid velocity, shows that the superfiuid velocity has a peak at the dip, and is directed 
inward, as we expect from the continuity equation. Thus, qualitatively, the onset of collapse 
happens on a time scale an order of magnitude greater than that for the establishment of 
the initial wave function, and it is initiated by a sudden inward rush of particles from a 
finite distance from the center. This signals the formation of a "black hole", as suggested 
by Ueda and Huang |]TB[. 
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After the onset of collapse, the density in the black hole increases rapidly, as indicated 
in Fig. §(c). The increase is fueled by particles drawn from throughout the outer region, 
but the density outside does not decrease uniformly. Instead, there are ripples of small 
implosions This is illustrated in the 3D plot of Fig. with r\ip\ plotted above the r-t plane. 
The black hole formation is indicated by the sudden rise of rl?/'! near the origin. The terrain 
is smooth before the rise, but very rough after that. The uneven bumps in the terrain are 
of course deterministic, but they are very sensitive to variations in the initial wave function, 
and in this sense they are random (see Fig. The fact that the height of the plateau seems 
to remain constant in time is an artifact of the finite spatial grid size. The height appears 
to increase without bound in the continuum limit, as illustrated in Fig. ^. 

To determine Nc more precisely, we plot the occurrence time of the first local minimum 
in the wave function within distance r = 0.05 from the origin, as a function of G = iV/2180. 
This is done in Fig. |^, for a Gaussian initial wave function with a = 1.0, and for an optimal 
Gaussian initial wave function with a = 1.79. The results are as follows: 

^ _ f 1145 (Gaussian initial state with a = 1.00) 
^ \ 1251 (Gaussian initial state with a = 1.79) 

which conforms to our expectation that the numbers are not very sensitive to initial condi- 
tions, because the time for the onset of collapse is long compared to the formation time of 
a quasi-stationary wave function. 

In Fig. 1^ we show the dependence of some of our results on the spatial grid size. The 
time for the onset of collapse appears to tend to a finite limit when the grid size approaches 
zero, but the height of a peak in the collapsed region seems to grow without bound, faster 
than (grid size)~^. This is consistent with the expectation that density fluctuations become 
divergent in the continuum limit. 



V. GROWTH-COLLAPSE CYCLES 



We have ignored gain-loss mechanisms so far, because we wanted to understand various 
effects one at a time. Now we shall take them into account, and make the problem more 
realistic. 

In the experimental situation, the condensate in the trap can exchange atoms with an 
uncondensed "thermal cloud" , which contains far more atoms than the condensate. Equilib- 
rium is reached when the chemical potentials equalize. The equilibrium fraction of conden- 
sate to "cloud" atoms is determined by the temperature. The kinetic equations governing 
the approach to equilibrium are quite intricate, and a subject of ongoing research ||21]]. We 
shall use a phenomenological gain equation based on a fit to data [E^]: 



dN 



coN 



N 



N, 



eq 



0.4 



(26) 



where A^eq is the equilibrium number of atoms in the condensate, with 



CO = /7 = /v^C(3/2)7ei(15Ge,)°-^T^ 

Kb J- 



(27) 
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where 7ei is the elastic scattering rate, which is approximately equal to 1 s~^ for the conditions 
of the experiment of Ref. [^] and C(3/2) ~ 2.612. / is a phenomenological factor to reconcile 
the above theoretical prediction for 7 with the experimental results for the repulsive case 



2^ . In our runs we used / = 5.75 and T = 300 nK, while Gen was taken to be equal to the 



critical value of G, i.e. about 0.574. Assuming N/N^q <C 1, we ignored the second term in 



the bracket. The main loss mechanisms are two- and three-body collisions p3| described by 
the following equation: 

dn s 2 

— = -cin - C2n 
at 

ci = 2.6 X 10"^^ cmVs 

C2 = 1.2 X 10"^^ cmVs (28) 

where n is the particle density in the condensate. Adding the gain and loss terms, we have 
an extended NLSE, which reads, in units such that = 2uj^^ = 1, 



du r 2 2G|mP\ / g2|m|4 gimP 



where 



70 = 2.6 X 10"^ 

71 = 8.3 X 10"^ 

72 = 7.2 X 10"^ (30) 

In Ref. 1 18], the value of 71 is taken to be of order 10""^, which differs considerably from 
ours. 

With the new equation, we find that the initial conditions become quite irrelevant. The 
number of particles quickly settle into a steady-state oscillation whose peak value defines 
the critical particle number in the new context. This is shown in Fig. 9, with: 

Nc ^ 1260 (31) 

The mean number of atoms is approximately 735, and the period of the oscillation is 

917 

r « = 1.01 s (32) 

The long time scale of the oscillations explains why they are insensitive to initial conditions. 
We have extended the calculations to the longest evolution time practicable, and found that, 
with a constant supply of particles, N oscillates in a steady state for at least 1.44x10^ trap 
cycles, or 100 seconds. These results are comparable with those obtained by Hulet et al. 
p!7| , except that, in contrast to our case, collapses to zero in that calculation. 
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FIGURES 



Condensate 
wave function 



Effective potential 




r > 

FIG. 1. Transfer a local value of the wave function from the left panel to the horizontal axis on 
the right. If it lands to the left of the potential barrier, the local system is classically stable, but 
can decay via tunneling, as indicated by the path A ^ B. If it lands to the right, the local system 
will rapidly collapse towards a state of infinite density. 
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FIG. 2. Time development of the condensate wave function. The time steps are given in units 
of 2uj^^, where lo is the circular frequency of the external harmonic trap. The number of particles 
is A'^ = 1142 (i.e. G = A^/2180 = 0.524), which is the highest particle number for which a Gaussian 
initial wave function with a = 1.0 is stable, (a) An initially uniform wave function tends to 
Gaussian form, except for local fluctuations, after a few time steps, (b) For an initial unperturbed 
Gaussian wave function with width do, the width narrows to a stationary value after a few time 
steps, (c) Starting from a Gaussian with width (io/\/1.48 gives a sensibly stationary wave function. 
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FIG. 3. Determining the best width d = do/s/a of the initial Gaussian wave function ipo{r), for 
N = 1248 {G = 0.5725), by minimizing the ratio of the standard deviation of the peak location to 
the mean peak location of |u(r, The minimum occurs at a = 1.79. 
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FIG. 4. Time development of unstable case, with A'" = 1380. (a) The wave function begins to 
narrow after t ~ 0.2. (b) The onset of collapse is signaled by the first occurrence of a dip (local 
minimum) in the wave function within a small distance from the origin, which indicates that a 
shell of particles begins to implode. The dip first appears at t 0.7 and moves quickly towards 
r = 0. After the onset, the collapse continues by drawing in particles from all distances, in ripples 
of implosive flow. In this figure the first 10 grid points (r ^ 0.05) were cut off for clarity, (c) The 
first local minimum of the wave function is now at less than r = 0.05. The wave function becomes 
more ragged as time goes on. Because of the finiteness of the grid the density at the center remains 
finite. 
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FIG. 5. Space-time behavior of the wave function of Fig. |^ in a 3D plot. The black hole formation 
is indicated by a sudden emergence of the high plateau. After that time, the terrain in the plain 
below the plateau becomes rough, indicating the ripples of implosive flow. As the spatial grid size 
tends to zero, the height of the plateau seems to increase without limit. (See Fig. ^.) 
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FIG. 6. Sensitivity of ripples to the initial conditions. We show the amplitude of the wave 
function just before and after the onset of the collapse for two nearly identical initial wave functions, 
Gaussians with a = 1 and a = 1.001. Notice that, before the collapse, the two wave functions 
practically coincide. However, once the collapse occurs, they diverge very quickly. 
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G=N/2180 

FIG. 7. Time for the onset of collapse vs Ar/2180. The two curves correspond to a Gaussian 

initial wave function with a = 1.0 and to a Gaussian initial wave function with a = 1.79. The 
critical number Nc is determined by the asymptote at which the time goes to infinity. We get 
Nc = 1145 (Gc = 00.525) when a = 1.0 and = 1251 (Gc = 0.574) when a = 1.79. The results 
for the critical number are not very sensitive to the initial wave function. 
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FIG. 8. Sensitivity of results on the spatial grid size, measured in units of cIq. The initial wave 
function is Gaussian with a = 1.0. Circles represent the time for the onset of collapse, which tends 
to a finite value in the continuum limit. Asterisks indicate the height of the peak of the square 
root of the density in the collapsed region. It appears to increase without limit as the grid size 
tends to zero. With the log-log axes used, the best linear fit has a slope of slightly less than — 1, 
which implies that the height of the peak increases faster than (grid size)~^ as the grid size goes 
to zero. 
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FIG. 9. Steady-state oscillation of the number of particles, when the condensate is fed by a 
"thermal cloud" and suffers loss through two- and three-body collisions. The peaks give Nc ~ 1260. 
The period of 1.01 s is much greater than the trap period of 7 ms. For this reason, initial conditions 
are irrelevant. 



18 



